OVERVIEW
Image credits: I-Hsuan Lin, National Yang-Ming University
OUR WORKFLOW
| PASSAGGIO | SOFTWARE | FORMATO FILE INPUT | FORMATO FILE OUTPUT | FILE SUPPLEMENTARI |
|---|---|---|---|---|
| Controllo qualità delle reads | FastQC | .fastq (.fq; compresso: .fq.gz) | .fastq (.fq; compresso: .fq.gz) | |
| Rimozione adattatori e pulizia basata su qualità | Trimmomatic | .fastq (.fq; compresso: .fq.gz) | .fastq (.fq; compresso: .fq.gz) | File fasta con sequenze adattatori |
| Mappaggio delle reads | STAR | .fastq (.fq; compresso: .fq.gz) | .sam / .bam (compresso) | File fasta del genoma e file GTF dell’annotazione |
| Visualizzazione reads mappate | IGV (Integrative Genomics Viewer) | .bam (indicizzato .bai) | Nessun output | File del genoma e annotazione (già disponibile su IGV) |
| Controllo qualità allineamento (opzionale) | Qualimap | .bam | Output tabulare | |
| Conta delle reads | featureCounts | .sam / .bam | .txt (testo) | File GTF dell’annotazione |
| Analisi geni differenzialmente espressi | DESeq2 / edgeR (pacchetti R, Bioconductor) | .txt (testo) | .txt / .csv (qualsiasi formato tabulare) | Tabella con disegno sperimentale |
| Data mining, Gene Ontology Enrichment Analysis (GOEA) | biomaRt, clusterprofiler | oggetti R derivati da analisi DE | grafici |
DAY 1 - 22 maggio 2018
RAW READS MANIPULATION
Accendere il computer e aprire il terminale (ctrl + alt + T)
Spostarsi nella cartella Scaricati, create una nuova cartella chiamandola tc
cd Scaricati/
mkdir tcEntrare nella cartella appena creata e valutarne il contenuto della cartella digitare il seguente comando
cd tc/
ls -lrthLa cartella dovrebbe essere vuota. Come prima cosa dobbiamo recuperare i file grezzi (FASTQ) del nostro esperimento. Per recuperare i file ci connettiamo al nostro server di laboratorio e con il protocollo di trasferimento ssh copiamo in locale (la cartella dove siamo in questo momento) i 2 file corrispondenti alle reads forward e alle reads reverse. Lavoreremo con un solo campione, ma come sapete l’esperimento consiste di 3 repliche biologiche per i non trattati e 3 repliche biologiche per quelli trattati. I comandi che qui eseguiremo possono essere ripetuti su più file contemporaneamente tramite appositi comandi come i cicli for
# prima scarichiamo il primo file (attenzione al punto in fondo al comando, indica current directory)
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/Sample_Untreated_1_1_subset.fastq.gz .
# inserire la password quando ce la chiede ()
# controlliamo che sia stato scaricato
ls -lrth
# ora scarichiamo il secondo file
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/Sample_Untreated_1_2_subset.fastq.gz .
# inserire nuovamente la password
# controlliamo
ls -lrthCome potete osservare la dimensione di questi file è molto bassa rispetto alla dimensione tipica di un file che risulta da un esperimento di sequenziamento come lo abbiamo descritto. Questo perché i file originali sono stati copiati e ridotti in dimensioni con lo scopo di riuscire a manipolarli e lavorarci nei tempi della lezione, principalmente per il fatto che le capacità dei computer che stiamo usando potrebbero essere limitate e quindi non in grado di eleaborare questi file in tempi brevi
Iniziamo con esplorare questi file
# "apriamo" il file
less -S Sample_Untreated_1_1_subset.fastq.gz
# contiamo il numero di reads
zcat Sample_Untreated_1_1_subset.fastq.gz | echo $((`wc -l`/4))
zcat Sample_Untreated_1_2_subset.fastq.gz | echo $((`wc -l`/4))
# il numero di reads F e R coincide?
FASTQ file
QUALITY CONTROL
Per effettuare il controllo di qualità sulle nostre reads utilizzeremo il programma FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Questo software (scritto in Java) prende come input un numero a piacimento di reads e calcola varie statistiche che possono dare un’idea della qualità dei nostri dati
# creiamo una nuova cartella dove scaricheremo il programma e salveremo i risultati
mkdir fastqc
cd fastqc/
# scarichiamo il programma con il comando wget
wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.7.zip
# scompattiamo la cartella
unzip fastqc_v0.11.7.zip
# rimuoviamo cartella compressa
rm fastqc_v0.11.7.zip
# vediamo cosa c'è dentro la cartella senza spostarci dentro (percorso relativo alla nostra posizione)
ls -lrth FastQC/
# il file che ci interessa, cioè il programma che andremo a richiamare è "fastqc"
# dobbiamo prima renderlo eseguibile
chmod 755 FastQC/fastqc
# controlliamo che funzioni
./FastQC/fastqc -hA questo punto possiamo lanciare il programma dando come input i due file con le reads F e R. Il comando è molto semplice, basta specificare i nomi dei file che vogliamo controllare, e la cartella dove vogliamo che vengano salvati i report coi risultati (la cartella deve già esistere quindi prima la creiamo)
# creiamo cartella per salvare i risultati
mkdir results
# ci spostiamo una cartella su, dove stanno i nostri file
cd ..
# eseguiamo il programma (dobbiamo specificare il percorso dell'eseguibile rispetto a dove siamo) e specifichiamo inoltre il percorso della cartella di output
./fastqc/FastQC/fastqc Sample_Untreated_1_1_subset.fastq.gz Sample_Untreated_1_2_subset.fastq.gz -o fastqc/results/Per ogni file di input vengono creati 2 file di output, un file html e una cartella compressa. Andiamo ad esplorare il report in formato html aprendolo con Firefox
FastQC report
A questo punto commentiamo insieme i risultati del controllo qualità confrontandoci con qualità di riferimento sul sito di FastQC e con la spiegazione dei moduli di analisi (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3%20Analysis%20Modules/)
Phred quality
Dopo aver deciso insieme quali parametri utilizzare per correggere la qualità delle nostre reads ed aver identificato gli adattatori da rimuovere, andiamo ad effettuare la rimozione degli adattatori e il quality trimming
ADAPTERS REMOVAL AND QUALITY TRIMMING
Per effettuare queste due operazioni useremo un unico programma che si chiama Trimmomatic (http://www.usadellab.org/cms/?page=trimmomatic). Il suo utilizzo è molto simile a quello di FastQC (anche Trimmomatic è scritto in Java) e necessita come input le sequenze grezze. Per ogni file di input (nel nostro caso 2 file) verranno creati 2 file di output, uno indicato come paired, l’altro come unpaired; è necessario indicare il nome di ciascuno dei file di output. Le opzioni di Trimmomatic sono:
PEindica che lavoriamo su file paired-end-trimlogindica che vogliamo salvare un file di log-threadsindica il numero di processori che vogliamo utilizzare per lanciare il programmaILLUMINACLIP:ALL_ADAPTERS.fa:2:30:10è la parte che indica la rimozione degli adattatori. Il file ALL_ADAPTERS.fa è il file in formato fasta che contiene la sequenza degli adattatori maggiormente utilizzatiLEADING:16indica la rimozione di N basi in testa alla sequenza, sotto la qualità indicataTRAILING:16indica la rimozione di N basi in coda alla sequenza, sotto la qualità indicataSLIDINGWINDOW:4:25effettua un trimming basato su una “finestra” di basi, rimuovendo le basi nel momento in cui la qualità media per quella finestra scende sotto la qualità indicataMINLEN:36indica la rimozione delle reads che sono più corte della lunghezza indicata
# ci spostiamo nella cartella tc, creiamo nuova cartella e scarichiamo il programma
# usiamo percorso assoluto da qualsiasi punto del computer per spostarci dentro una qualsiasi cartella
cd /home/usr/Scaricati/tc # usr è specifico per il vostro utente (il mio è dlfptr70)
mkdir Trimmomatic
wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.36.zip -O Trimmomatic/Trimmomatic-0.36.zip # con l'opzione -O specifichiamo dove vogliamo scaricare il file e con che nome
# decomprimiamo senza cambiare cartella e rimuoviamo l'archivio
unzip Trimmomatic/Trimmomatic-0.36.zip -d Trimmomatic/ && rm Trimmomatic/Trimmomatic-0.36.zip
# controlliamo il funzionamento del programma
java -jar Trimmomatic/Trimmomatic-0.36/trimmomatic-0.36.jar -h # -h è l'opzione per helpA questo punto possiamo lanciare Trimmomatic sui nostri file; anche se dal controllo qualità non abbiamo identificato sequenze di adattatori come contaminanti, un classico approccio potrebbe essere quello di dare come input di file degli adattatori il file fornito dal programma che quindi rimuoverà tutte le sequenze corrispondenti a quelle che trova
# creiamo file unendo in uno tutti gli adattatori
cat Trimmomatic/Trimmomatic-0.36/adapters/*.fa > Trimmomatic/Trimmomatic-0.36/adapters/ALL_ADAPTERS.fa
# cat è un comando che unisce file "row-wise"; l'asterisco (*.fa) indica di selezionare tutti i file che hanno estensione .fa
java -jar Trimmomatic/Trimmomatic-0.36/trimmomatic-0.36.jar PE -threads 10 -trimlog \
Trimmomatic/log_trimmomatic \
Sample_Untreated_1_1_subset.fastq.gz Sample_Untreated_1_2_subset.fastq.gz \
Trimmomatic/Sample_Untreated_1_1_subset_trim_paired.fastq.gz Trimmomatic/Sample_Untreated_1_1_subset_trim_unpaired.fastq.gz \
Trimmomatic/Sample_Untreated_1_2_subset_trim_paired.fastq.gz Trimmomatic/Sample_Untreated_1_2_subset_trim_unpaired.fastq.gz \
ILLUMINACLIP:Trimmomatic/Trimmomatic-0.36/adapters/ALL_ADAPTERS.fa:2:30:10 \
LEADING:16 TRAILING:16 SLIDINGWINDOW:4:25 MINLEN:36 &> Trimmomatic/log_trimmomatic_screen
# con il simbolo &> andiamo a salvare su un file anche l'output che normalmente sarebbe stampato sullo schermoDovrebbe terminare dopo qualche secondo senza aver stampato a schermo alcuna informazione. Andiamo a vedere se ha creato i file ed esaminiamo i file di log (log_trimmomatic e log_trimmomatic_screen); il primo contiene le seguenti info per ciascuna read:
- the read name
- the surviving sequence length
- the location of the first surviving base, aka. the amount trimmed from the start
- the location of the last surviving base in the original read
- the amount trimmed from the end
Il secondo contiene le statistiche generali
# guardiamo prima le statistiche generali
less -S Trimmomatic/log_trimmomatic_screen
# poi l'altro file
less -S Trimmomatic/log_trimmomaticA questo punto ripetiamo il controllo di qualità, questa volta sulle reads trimmate, solo sui file paired
# creiamo cartella dove salveremo nuovo controllo qualità
mkdir fastqc/results_after_trimming
# lanciamo FastQC sulle sequenze trimmate
./fastqc/FastQC/fastqc Trimmomatic/Sample_Untreated_1_1_subset_trim_paired.fastq.gz \
Trimmomatic/Sample_Untreated_1_2_subset_trim_paired.fastq.gz -o fastqc/results_after_trimming/Andiamo ad esplorare i report e vedere se le statistiche sulla qualità sono cambiate
MAPPING
Ora i nostri file sono pronti per essere allineati sul genoma di riferimento (Arabidopsis). Per quasi tutti i programmi che effettuano allineamento il primo passaggio consiste della creazione dell’indice del genoma. Perché viene creato l’indice del genoma: citazione presa da qualcuno ma non ricordo chi (modificata) “Se entro in una biblioteca con la frase “La genetica è il ramo della biologia che si occupa del materiale ereditario” e voglio trovare il libro in cui è scritta, farò molta fatica a trovarlo, a meno che all’ingresso della biblioteca non si trovi un catalogo che mi indirizza almeno verso la sezione “Scienza””. L’indice quindi consiste in una serie di file che riassumono le principali caratteristiche del genoma di un organismo. Per effettuare l’allineamento useremo un programma chiamato STAR (https://github.com/alexdobin/STAR). Andiamo a scaricarlo insieme ai file aggiornati del genoma e del trascrittoma di Arabidopsis
# siamo nella cartella tc
mkdir STAR
# scarichiamo e salviamo nella cartella appena creata
wget https://github.com/alexdobin/STAR/archive/2.6.0c.tar.gz -O STAR/2.6.0c.tar.gz
# decomprimiamo e rimuoviamo archivio
tar -xf STAR/2.6.0c.tar.gz -C STAR/ && rm STAR/2.6.0c.tar.gz
# controlliamo il funzionamento
./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STAR -hPer creare l’indice del genoma, abbiamo bisogno del genoma. Andiamo a scaricarlo dal sito di Ensembl Plants
# genoma
mkdir genome
wget ftp://ftp.ensemblgenomes.org/pub/release-39/plants/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz -O genome/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
# trascrittoma
wget ftp://ftp.ensemblgenomes.org/pub/release-39/plants/gtf/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.39.gtf.gz -O genome/Arabidopsis_thaliana.TAIR10.39.gtf.gz
# i file devono essere decompressi
gunzip genome/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
gunzip genome/Arabidopsis_thaliana.TAIR10.39.gtf.gzA questo punto possiamo creare l’indice con il comando di STAR --runMode genomeGenerate. Vediamo nel dettaglio le altre opzioni:
./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STARchiamiamo il programma--runThreadNnumero di threads con cui vogliamo lanciare il programma--runModeindica che vogliamo creare l’indice--genomeDirindica la cartella dove vogliamo che venga salvato l’output--genomeFastaFilescartella dove è presente il file del genoma--sjdbGTFfileindichiamo che vogliamo usare anche il trascrittoma per aiutarlo nella creazione dell’indice--sjdbOverhangspecifica la lunghezza delle reads meno 1. Serve per creare il database delle giunzioni di splicing
Lanciamo il comando con le opzioni e i percorsi adattati al nostro caso
# creiamo cartella per l'indice all'interno della cartella del genoma
mkdir -p genome/index
# lanciamo la creazione dell'indice
./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STAR \
--runThreadN 10 \
--runMode genomeGenerate \
--genomeDir genome/index \
--genomeFastaFiles genome/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa \
--sjdbGTFfile genome/Arabidopsis_thaliana.TAIR10.39.gtf \
--sjdbOverhang 100Con 10 core dovrebbe impiegare circa 3 minuti. Andiamo a vedere se ha creato i file nella cartella index
# spostiamo il file Log.out dove dovrebbe stare
mv Log.out genome/index/
# guardiamo i file appena creati
ls -lrht genome/index
# raccolgono diverse statistiche sul genoma di Arabidopsis, comprese info sul db di splicejunctionsPossiamo ora andare ad effettuare l’allineamento. Ovviamente per l’allineamento useremo le reads trimmate e pulite (le paired). Vediamo nel dettaglio i comandi:
./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STARchiamiamo il programma--runThreadNnumero di threads con cui vogliamo lanciare il programma--genomeDirdove sta la cartella con l’indice--readFilesInqui scriviamo il percorso e il nome dei campioni da mappare--readFilesCommandqual è il formato di comrpessione dei file da mappare--outFileNamePrefixpercorso e prefisso che vogliamo assegnare ai file di output--alignIntronMaxlunghezza massima (stimata) degli introni del nostro organismo--outSAMtypetipo es estensione del file di output
Attenzione che la maggior parte dei programmi di allineamento sono creati con impostazioni di default basate sul genoma umano –> aggiustare i parametri sulla base del proprio organismo! Andiamo a lanciare il comando
# creiamo cartella dove salveremo file di allineamento
mkdir -p STAR/aligned
# lanciamo l'allineamento
./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STAR \
--runThreadN 10 \
--genomeDir genome/index/ \
--readFilesIn Trimmomatic/Sample_Untreated_1_1_subset_trim_paired.fastq.gz \
Trimmomatic/Sample_Untreated_1_2_subset_trim_paired.fastq.gz \
--readFilesCommand zcat \
--outFileNamePrefix STAR/aligned/SU_1_ \
--alignIntronMax 2000 \
--outSAMtype BAM SortedByCoordinateSe è tutto a posto l’allineamento dovrebbe impiegare circa 10 secondi. Come vedete la lunghezza massima media per gli introni di Arabidopsis è di circa 2000 bp, e gli ho detto che voglio che il file di output sia in formato .bam (formato allineamento binario). Il file .bam è in un formato che non può essere aperto con il terminale, ma che andremo poi a visualizzare con un apposito software. Per completezza vi inserisco un immagine dello stesso file di allineamento ma in un formato visualizzabile e consultabile sul terminale. Il formato di questo file è .sam (Sequence Alignment Format), formato non molto utilizzato perché non essendo compresso occupa molto spazio
SAM file
La prima riga descrive il numero della versione del file
Dalla seconda alla ottava sono riportate le lunghezze in paia di basi dei cromosomi di Arabidopsis
La nona e la decima riga indicano il programma utilizzato per l’allineamento e il comando digitato
Dalla undicesima riga ogni riga coincide con una read, le colonne sono descritte di seguito
La colonna più interessante è forse la sesta (CIGAR) dove vengono specificati nel dettaglio gli allineamenti di ciascuna reads con il reference
SAM file description
Andiamo ora a vedere i nostri file creati
# il file che ci interessa è chiamato SU_1_Aligned.sortedByCoord.out.bam
ls -rht STAR/aligned/
# leggiamo qualche statistica in questo altro file
less -S STAR/aligned/SU_1_Log.final.outQuanto vi risulta la percentuale di Uniquely mapped reads?
ALIGNMENT VISUALIZATION
Dopo aver effettuato l’allineamento può essere utile effettuare un controllo visivo sull reads allineate. A tal fine si utilizza spesso un software chiamato IGV (http://software.broadinstitute.org/software/igv/) che permette tra le altre cose di visualizzare i file dell’allineamento dopo aver caricato un genoma di riferimento
Integrative Genomics Viewer
Andiamo a scaricare il programma dal sito http://software.broadinstitute.org/software/igv/download cliccando in basso su “Launch with 2 GB”. Attendiamo lo scaricamento e apriamo il programma con doppio clic. Una volta aperto in alto a sinistra selezionare il genoma di A. thaliana TAIR10. A questo punto possiamo caricare i risultati del nostro allineamento, non prima però di aver creato un indice per il file .bam. Un altro indice direte voi, sì questa volta però non eseguiremo nessun passaggio ma vi farò prelevare il file dell’indice già preparato dal nostro server. Per completezza vi mostro il comando utilizzato per creare l’indice con Samtools
samtools index -b SU_1_Aligned.sortedByCoord.out.bamSu questi computer non possiamo installare il programma che si utilizza, tra le altre cose, anche per la creazione degli indici per i file di allineamento. Il programma si chiama Samtools (http://www.htslib.org/) ed è molto utilizzato in bioinformatica perché svolge numerose funzioni. Andiamo a copiare nella cartella dove sta il file .bam appena creato da noi, il file dell’indice, che permette che il file .bam venga visualizzato correttamente in IGV
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/SU_1_Aligned.sortedByCoord.out.bam.bai STAR/aligned/
# il file che vedete deve avere lo stesso nome del file di allineamento, ma ha un'estensione diversa, .baiA questo punto potete tornare su IGV e cliccare su “File -> Load from file” e selezionare dalla cartella aligned il file SU_1_Aligned.sortedByCoord.out.bam. Selezioniamo Chr1 e zoomiamo sino a che la rotella inizia a girare, a quel punto possiamo visualizzare le nostre reads allineate ai geni del genoma di riferimento
Non sarà così facile trovare un po’ di reads!
Nel caso non funzionasse, copiare dal server alla vostra cartella, sia il file .bam che il .bai (potrebbe essere che il bai debba essere creato dallo stesso bam per funzionare la visualizzazione)
# bam
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/SU_1_Aligned.sortedByCoord.out.bam STAR/aligned/
# bai
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/SU_1_Aligned.sortedByCoord.out.bam.bai STAR/aligned/La lunghezza degli introni è stata rispettata?
DAY 2 - 29 maggio 2018
READS COUNTING
Una volta controllato che l’allineamento è andato a buon fine e siamo soddisfatti di come il programma di mapping ha lavorato, possiamo spostarci nel prossimo step verso l’analisi dei geni, ovvero la conta delle reads. Per conta delle reads si intende il contare quante reads sono state assegnate (hanno “mappato”) a un determinato gene (o trascritto, o esone, le cosiddette genomic features) in seguito al mapping. Per questo passaggio utilizzeremo il software featureCounts (http://subread.sourceforge.net/) del pacchetto Subread. Scarichiamo e testiamo il programma
# nella cartella tc/ creiamo cartella
mkdir featureCounts
# scarichiamo dal server perché più comodo
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/featureCounts featureCounts/
# testiamo funzionamento
./featureCounts/featureCounts -vVediamo le opzioni e i parametri di featureCounts:
featureCountslanciamo il programma-aforniamo il file di annotazione del trascrittoma-opercorso e nome file di output-pindichiamo che lavoriamo su reads paired e che verranno contati i frammenti e non le reads-tindica quale feature contare-gindica che vogliamo raggrppare tutte le feature in questa meta-feature-Tquanti threads vogliamo usare per correre il programmain.bamnome del file di allineamento che vogliamo come input
A questo punto andiamo a lanciare il programma, ricordandoci che nel nostro caso vogliamo contare le feature del file di allineamento .bam
# creiamo cartella per output
mkdir -p featureCounts/output
# lanciamo la conta
./featureCounts/featureCounts \
-a genome/Arabidopsis_thaliana.TAIR10.39.gtf \
-o featureCounts/output/SU_1_counts.txt \
-p \
-t exon \
-g gene_id \
-T 10 \
STAR/aligned/SU_1_Aligned.sortedByCoord.out.bamDovrebbe terminare in qualche secondo. Andiamo ad esplorare i file di output
ls -lrht featureCounts/output/
less -S featureCounts/output/SU_1_counts.txt.summary
less -S featureCounts/output/SU_1_counts.txt
# vediamo solo le info che ci interessano
cut -f 1,7 featureCounts/output/SU_1_counts.txt | sed 1d | less -SLe colonne che ci interessano sono la prima e l’ultima (7). Perché per alcuni geni abbiamo più di un valore per le colonne 2, 3, 4 e 5? I passaggi che abbiamo visto finora in una situazione normale andrebbero ripetuti per tutti i campioni che abbiamo fatto sequenziare, sino ad ottenere una matrice di conte, come l’ultimo file che abbiamo visto, per ciascun campione. Da questa matrice di conte (counts matrix) creeremo un file unico che sarà l’input per la fase successiva. Visto che nella nostra esercitazione abbiamo elaborato solo un file, vi fornirò io la counts matrix completa (6 campioni, 3 controlli e 3 trattati)
# scarichiamo dal server il file con le conte per i 3 controlli e i 3 trattati
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/countDF_unvs3h featureCounts/
# vediamo il file
less -S featureCounts/countDF_unvs3hINSTALLING R AND RSTUDIO
A questo punto abbiamo concluso la prima fase del corso che abbiamo svolto usando solamente il terminale. Per la seconda fase utilizzeremo invece il software R e il suo ambiente di sviluppo integrato (IDE) chiamato RStudio. Per prima cosa vi invito a scaricare e installare R ed RStudio seguendo le istruzioni di seguito a seconda del vostro SO:
Windows: cliccare su http://cran.mirror.garr.it/mirrors/CRAN/bin/windows/base/R-3.5.0-win.exe, attendere lo scaricamento e installare il file con estensione
.exeappena scaricato, lasciare le impostazioni di default. Per installare RStudio andare alla pagina https://www.rstudio.com/products/rstudio/download/#download e alla voce “Installers for Supported Platforms” cliccare su “RStudio 1.1.453 - Windows Vista/7/8/10”, attendere lo scaricamento e installare il file con estensione.exeappena scaricato, lasciare le impostazioni di default.MacOS: cliccare su http://cran.mirror.garr.it/mirrors/CRAN/bin/macosx/R-3.5.0.pkg, attendere lo scaricamento e installare il file con estensione
.pkgappena scaricato, lasciare le impostazioni di default. Per installare RStudio andare alla pagina https://www.rstudio.com/products/rstudio/download/#download e alla voce “Installers for Supported Platforms” cliccare su “RStudio 1.1.453 - Mac OS X 10.6+ (64-bit)”, attendere lo scaricamento e installare il file con estensione.dmgappena scaricato, lasciare le impostazioni di default.Linux: informazioni disponibili all’indirizzo https://medium.com/@GalarnykMichael/install-r-and-rstudio-on-ubuntu-12-04-14-04-16-04-b6b3107f7779
Per qualsiasi problema scrivetemi all’indirizzo delfino.pietro@gmail.com o possiamo vedere domani l’installazione in aula, ma sarebbe meglio arrivare già con entrambi i programmi installati e alcuni pacchetti che elenco di seguito, con i rispettivi comandi per l’installazione. Aprire RStudio e nella console digitare i seguenti comandi con sfondo grigio, in ordine, uno per uno:
- Bioconductor
source("https://bioconductor.org/biocLite.R")biocLite()
- tidyverse
install.packages("tidyverse")
- DESeq2
source("https://bioconductor.org/biocLite.R")biocLite("DESeq2")
- edgeR
source("https://bioconductor.org/biocLite.R")biocLite("edgeR")
- biomaRt
source("https://bioconductor.org/biocLite.R")biocLite("biomaRt")
- clusterProfiler
source("https://bioconductor.org/biocLite.R")biocLite("clusterProfiler")
- pheatmap
install.packages("pheatmap")
- RColorBrewer
install.packages("RColorBrewer")
- Genome wide annotation of Arabidopsis
source("https://bioconductor.org/biocLite.R")biocLite("org.At.tair.db")
DIFFERENTIAL GENE EXPRESSION ANALYSIS
DESEQ2
Torniamo alla lezione. Apriamo RStudio e ci spostiamo nella cartella dove abbiamo il file completo con le conte
setwd("/home/usr/Scaricati/tc/featureCounts")
# "usr" è variabile a seconda del vostro nome utente, usate il tab per completare il percorso# carichiamo la libreria
library(DESeq2)
# carichiamo il file
cts <- as.matrix(read.table("featureCounts/countDF_unvs3h", sep="\t",
header=T, col.names = c("ID","un1","un2","un3","NO3_1","NO3_2","NO3_3"),
row.names = "ID"))# non eseguire questo comando, è solo per me
cts <- as.matrix(read.table("Scaricati/tc/featureCounts/countDF_unvs3h", sep="\t",
header=T, col.names = c("ID","un1","un2","un3","NO3_1","NO3_2","NO3_3"),
row.names = "ID"))
# vediamo cosa abbiamo caricato
head(cts)## un1 un2 un3 NO3_1 NO3_2 NO3_3
## AT1G01010 248 535 190 830 839 497
## AT1G01020 255 684 943 842 679 909
## AT1G01030 546 607 568 61 45 84
## AT1G01040 690 1970 1386 607 687 419
## AT1G01046 2 7 7 1 3 0
## AT1G01050 597 1252 2008 1015 1054 1061# creiamo oggetto per il disegno sperimentale
coldata <- data.frame(row.names = colnames(cts))
coldata$condition <- c(rep("untreated", 3), rep("treated", 3))
coldata$time <- c(rep("0", 3), rep("3", 3))
str(coldata)## 'data.frame': 6 obs. of 2 variables:
## $ condition: chr "untreated" "untreated" "untreated" "treated" ...
## $ time : chr "0" "0" "0" "3" ...# convertiamo le colonna in factor per poterle usare nella formula
coldata$condition <- as.factor(coldata$condition)
coldata$time <- as.factor(coldata$time)
str(coldata)## 'data.frame': 6 obs. of 2 variables:
## $ condition: Factor w/ 2 levels "treated","untreated": 2 2 2 1 1 1
## $ time : Factor w/ 2 levels "0","3": 1 1 1 2 2 2# costruiamo l'oggetto DESeq
dds <- DESeqDataSetFromMatrix(countData = cts,
colData = coldata,
design = ~ condition )
# ordiniamo i livelli per usare "untreated" come confronto
dds$condition <- relevel(dds$condition, ref = "untreated")
# vediamo info sull'oggetto creato
dds## class: DESeqDataSet
## dim: 33602 6
## metadata(1): version
## assays(1): counts
## rownames(33602): AT1G01010 AT1G01020 ... ATMG01400 ATMG01410
## rowData names(0):
## colnames(6): un1 un2 ... NO3_2 NO3_3
## colData names(2): condition time# rimuoviamo geni mai espressi
dds <- dds [rowSums(counts(dds)) > 0, ]
dds## class: DESeqDataSet
## dim: 24383 6
## metadata(1): version
## assays(1): counts
## rownames(24383): AT1G01010 AT1G01020 ... ATMG01400 ATMG01410
## rowData names(0):
## colnames(6): un1 un2 ... NO3_2 NO3_3
## colData names(2): condition time# quanti erano?Data exploration with DESeq2
library(tidyverse)
# trasformiamo i dati per la PCA; correggiamo per rendere dataset omoschedastico
# facciamo 2 trasformazioni proprie di DESeq2, rlog e vst
rld <- rlog(dds, blind = FALSE)
head(assay(rld),3)## un1 un2 un3 NO3_1 NO3_2 NO3_3
## AT1G01010 8.878895 8.722974 7.173810 9.572869 9.646944 8.948716
## AT1G01020 8.986913 9.101649 9.220015 9.659787 9.448712 9.779084
## AT1G01030 9.716037 8.717315 8.348288 6.198081 5.911402 6.591842vsd <- vst(dds, blind = FALSE)
head(assay(vsd),3)## un1 un2 un3 NO3_1 NO3_2 NO3_3
## AT1G01010 8.992480 8.825655 7.253611 9.743679 9.824669 9.067229
## AT1G01020 9.029901 9.155565 9.283881 9.763533 9.532699 9.894374
## AT1G01030 10.075139 8.994522 8.605132 6.542427 6.302154 6.885296# calcoliamo PCA sui dati trasformati con rlog usando la funzione di R,
# dist, che calcola distanze Euclidee sulla matrice trasposa rld
sampleDists <- dist(t(assay(rld)))
sampleDists## un1 un2 un3 NO3_1 NO3_2
## un2 58.79275
## un3 75.19080 64.52936
## NO3_1 199.18353 195.74731 200.74712
## NO3_2 197.94528 194.40729 198.29539 45.88588
## NO3_3 185.94000 186.16615 177.45544 93.93457 93.37848# plottiamo heatmap delle distanze specificando noi stessi la distanza che abbiamo
# appena calcolato con l'argomento "clustering distance"
sampleDistMatrix <- as.matrix( sampleDists )
rownames(sampleDistMatrix) <- rownames(coldata)
colnames(sampleDistMatrix) <- rownames(coldata)
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
clustering_distance_rows = sampleDists,
clustering_distance_cols = sampleDists,
col = colors
)
# plottiamo PCA
# custom PCA plot
pcaData <- plotPCA(rld, intgroup=c("condition", "time"), returnData=TRUE)
pcaData.vsd <- plotPCA(vsd, intgroup=c("condition", "time"), returnData=TRUE)
ggplot(pcaData, aes(PC1, PC2, color=condition, shape=time)) +
geom_text_repel(aes(label=pcaData$name)) +
geom_point(size=3) +
theme_bw() + ggtitle("Data transformed with rld")
ggplot(pcaData.vsd, aes(PC1, PC2, color=condition, shape=time)) +
geom_text_repel(aes(label=pcaData$name)) +
geom_point(size=3) +
theme_bw() + ggtitle("Data transformed with vsd")
Dopo alcune analisi esplorative possiamo eseguire il test della differenza di espressione tra i geni vero e proprio. Teniamo presente che DESeq2 lavora sulle conte grezze per effettuare il test, e non sui dati trasformati, quindi useremo l’oggetto dds e non le trasformazioni rld o vsd
# comando che effettua l'analisi vera e propria
dds <- DESeq(dds)
# salviamo oggetto con risultati specificando quali confronti
res <- results(dds, contrast=c("condition","treated","untreated"), alpha = 0.01)
summary(res)##
## out of 24383 with nonzero total read count
## adjusted p-value < 0.01
## LFC > 0 (up) : 4570, 19%
## LFC < 0 (down) : 4445, 18%
## outliers [1] : 11, 0.045%
## low counts [2] : 3782, 16%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results# we can order our results table by the smallest p-value
resOrdered <- res[order(res$pvalue),]
resOrdered## log2 fold change (MLE): condition treated vs untreated
## Wald test p-value: condition treated vs untreated
## DataFrame with 24383 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue
## <numeric> <numeric> <numeric> <numeric> <numeric>
## AT1G23740 5381.069 -5.994458 0.1754464 -34.16690 7.503470e-256
## AT1G22400 24333.995 7.492121 0.2250072 33.29724 4.232379e-243
## AT2G04040 9615.186 9.989993 0.3080733 32.42733 1.130964e-230
## AT2G21640 8584.189 6.425284 0.1986819 32.33956 1.945127e-229
## AT2G04050 8942.769 9.455977 0.3116858 30.33817 3.598939e-202
## ... ... ... ... ... ...
## AT4G25470 773.59979 0.1412803 1.313684 0.1075451 NA
## AT5G35935 35.60476 4.6814514 2.768811 1.6907802 NA
## AT5G37770 1252.46550 0.2171257 1.240675 0.1750061 NA
## AT5G44430 872.74506 -1.0961080 2.117555 -0.5176291 NA
## AT5G62920 737.38478 -1.2402995 1.353632 -0.9162750 NA
## padj
## <numeric>
## AT1G23740 1.544964e-251
## AT1G22400 4.357235e-239
## AT2G04040 7.762185e-227
## AT2G21640 1.001254e-225
## AT2G04050 1.482043e-198
## ... ...
## AT4G25470 NA
## AT5G35935 NA
## AT5G37770 NA
## AT5G44430 NA
## AT5G62920 NAVolcano plot
# Threshold on the adjusted p-value
alpha <- 0.01
# colors for the points
cols <- densCols(res$log2FoldChange, -log10(res$pvalue))
# plot
plot(res$log2FoldChange, -log10(res$padj), col=cols, panel.first=grid(),
main="Volcano plot", xlab="Effect size: log2(fold-change)", ylab="-log10(adjusted p-value)",
pch=20, cex=0.6)
abline(v=0)
abline(v=c(-2,2), col="brown")
abline(h=-log10(alpha), col="brown")